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We present the algorithmic details of the dynamical cluster 
approximation (DCA), with a quantum Monte Carlo (QMC) 
method used to solve the effective cluster problem. The DCA 
is a fully-causal approach which systematically restores non- 
local correlations to the dynamical mean field approximation 
(DMFA) while preserving the lattice symmetries. The DCA 
becomes exact for an infinite cluster size, while reducing to 
the DMFA for a cluster size of unity. We present a generaliza- 
tion of the Hirsch-Fye QMC algorithm for the solution of the 
embedded cluster problem. We use the two-dimensional Hub- 
bard model to illustrate the performance of the DCA tech- 
nique. At half-filling, we show that the DCA drives the spuri- 
ous finite-temperature antiferromagnetic transition found in 
the DMFA slowly towards zero temperature as the cluster 
size increases, in conformity with the Mermin- Wagner the- 
orem. Moreover, we find that there is a finite temperature 
metal to insulator transition which persists into the weak- 
coupling regime. This suggests that the magnetism of the 
model is Heisenberg like for all non-zero interactions. Away 
from half-filling, we find that the sign problem that arises in 
QMC simulations is significantly less severe in the context of 
DCA. Hence, we were able to obtain good statistics for small 
clusters. For these clusters, the DCA results show evidence 
of non- Fermi liquid behavior and superconductivity near half- 
filling. 

I. INTRODUCTION 

One of the most active subfields in condensed matter 
theory is the development of new algorithms to simu- 
late the many-body problem. This interest is motivated 
by various physical phenomena, including high tempera- 
ture superconductivity, magnetism, heavy fermions and 
the rich phenomenology occurring in quasi-one dimen- 
sional compounds. In the last few years, important 
progress has been made. Well-controlled results have 
been obtained by exact diagonalization and quantum 
Monte Carlo methods (QMC) However, these algo- 
rithms suffer from a common limitation in that the num- 
ber of degrees of freedom grows rapidly with the lattice 
size. As a consequence, the calculations are restricted 
to relatively small systems. In most cases, the limited 
size of the system prohibits the study of the low-energy 
physics of these models. 

Recently, another route to quantum simulations has 
been proposed. Following Metzner and VoUhardt and 
Miiller-Hartmann Q who showed that in the limit of infi- 
nite dimensions, the many-body problem becomes purely 
local, a mapping to a self-consistent Anderson impurity 



problem was performed The availability of many 

techniques to solve the Anderson impurity Hamiltonian 
has led to a dramatic burst of activity. However, when 
applied to systems in two or three dimensions this self- 
consistent approximation, referred to as the dynamical 
mean field approximation (DMFA) , displays some limita- 
tions. Due to its local nature, the DMFA neglects spatial 
fluctuations which are essential when the order parame- 
ter is non-local, or when short-ranged spin correlations 
are present. 

An acceptable theory which systematically incorpo- 
rates non-local corrections to the DMFA is needed. It 
must be able to account for fluctuations in the local 
environment in a self-consistent way, become exact in 
the limit of large cluster sizes, and recover the DMFA 
when the cluster size equals one. It must be easily im- 
plementable numerically and preserve the translational 
and point-group symmetries of the lattice. Finally, it 
should be fully causal so that the single-particle Green 
function and self energy are analytic in the upper half 
plane. There have been several attempts to formulate 
theories which satisfy these requirements, but all fail in 
some significant way 1^. 

In recent publications , the dynamical cluster ap- 
proximation (DCA) has been proposed as an extension 
to the DMFA which satisfies all these criteria. The DCA 
is built in close analogy with the DMFA. In the DCA, 
the lattice problem is mapped to a self-consistently em- 
bedded finite-sized cluster, instead of a single impurity 
as in the DMFA. The key idea of the DCA is to use the 
irreducible quantities (self energy, irreducible vertices) of 
the embedded cluster as an approximation for the corre- 
sponding lattice quantities. These irreducible quantities 
are then applied to construct the lattice reducible quan- 
tities such as the Green function or susceptibilities in the 
different channels. The cluster problem generated by the 
DCA may be solved by using a variety of techniques in- 
cluding the Quantum Monte Carlo (QMC) metho d 
the Fluctuation Exchange (FLEX) approximation [[11[ or 
the Non-Crossing Approximation (NCA) |p^ . 

The QMC method, and the Hirsch-Fye algorithm |l^] 
in particular, is the most reliable of these techniques. 
The Hirsch-Fye algorithm was originally designed for the 
treatment of few-impurity problems. Hence, it has been 
widely applied to the Kondo problem ||l^ and also to 
solve the impurity problem of the DMFA. For embedded 
cluster problems, this algorithm shows a mild sign prob- 
lem, compared to that encountered in previous finite- 
sized simulations, presumably due to the coupling to the 
host. Thus, we are able to perform simulations at sig- 
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nificantly lower temperatures than with other available 
techniques. However, in order to study a meaningful set 
of clusters of different sizes, it is necessary to use mas- 
sively parallel computers. 

Throughout this paper, we will use the two-dimensional 
Hubbard model on a simple square lattice as an example. 
Its Hamiltonian is 

H = - ^U-i{4^Cj„ + \i.c.) + e'^rii^ 

ij io- 

+ C/^KT-l/2)(nu~f/2) (1) 

i 

where is the matrix of hopping integrals, c[P is the 
annihilation (creation) operator for electrons on lattice 
site i with spin ct, and U the intraatomic repulsion. We 
will take /i = and vary the orbital energy e to fix 
the filling. The model has a long history and is still 
the subject of an intensive research in relation with the 
high-temperature superconductivity, the non-Fermi liq- 
uid phenomenon, the metal to insulator transition and 
magnetism in various physical systems dominated by 
strong correlations. Some short accounts |f gjlf j| ] of the 
DCA applied to this model have been recently published 
but without a full description of the details of the algo- 
rithm and numerical subtleties. It is the purpose of this 
paper to present the full account of the DCA-QMC tech- 
nique. A typical DCA algorithm using the QMC tech- 
nique as the cluster solver is made of three main blocks: 
the self-consistent loop, the analysis block and the an- 
alytical continuation block. The self-consistent loop is 
the most important of the three blocks; it is devoted to 
the mapping of the lattice to the cluster (coarse-graining) 
and to the solution of the cluster problem by the QMC 
method. In the analysis block, cluster Green functions 
obtained from the QMC method are transformed to lat- 
tice Green functions following the procedure described in 
section III . The last block is devoted to the computation 
of the lattice real frequency quantities from the analytical 
continuation of the corresponding QMC imaginary-time 
quantities by the maximum entropy method (MEM). ||l^] 
This paper is organized as follows. In the next section, 
we review the dynamical mean field approximation. In 
Sec. [II, we review the DCA formalism in which the lat- 



tice problem is mapped to a self-consistently embedded 
periodic cluster, and discuss the relationship between the 
cluster and the lattice. In this section we describe how 
different lattice Green functions can be obtained from the 



cluster quantities. In Sec. IV, we derive a modified form 
of the Hirsch-Fye QMC algorithm, which may be used to 
solve the effective cluster problem. We also discuss the 
conditioning and optimization of a variety of one and two- 
particle measurements. In Sec. ^ we discuss the DCA 



two techniques which has been discussed in earlier pub- 
lications [ p^Ul^Jl^ . At half-filling we discuss the occur- 
rence of antifcrromagnetism and the metal to insulator 
transition. Away from half-filling, we show the signature 
of a non- Fermi liquid behavior and superconductivity for 
small clusters for which the negative sign problem is mild 
so that good statistics can be obtained at low tempera- 
tures. Finally, in Sec. VII, we draw the conclusions on 
the present work and discuss future applications of the 
DCA to various physical problems. 



II. THE DYNAMICAL MEAN-FIELD 
APPROXIMATION 

The DCA algorithm is constructed in analogy with the 
DMFA. The DMFA is a local approximation which was 
used by various authors in perturbative calculations as 
a simplification of the k-summations which render the 
problem intractable 18|. But it was after the work 
of Metzner and VoUhardt ||^ and Miiller-Hartmann j|] 
who showed that this approximation becomes exact in 
the limit of infinite dimension that it received extensive 
attention during the last decade. In this approximation, 
one neglects the spatial dependence of the self-energy, 
retaining only its variation with time. Please see the 
reviews by Pruschke et al ^ and Georges et al |^ for a 
more extensive treatment. In this section, we will show 
that it is possible to re-interpret the DMFA as a course 
graining approximation, and then review its derivation. 

The DMFA consists of mapping the original lattice 
problem to a self-consistent impurity problem. As il- 
lustrated in Fig. |l| for a two-dimensional lattice, this is 
equivalent to averaging the Green functions used to cal- 
culate the irreducible diagrammatic insertions over the 
Brillouin zone. An important consequence of this aver- 
aging is that the self-energy and the irreducible vertices 
of the lattice are independent of the momentum. Hence, 
they are those of the impurity. 



algorithm. In Sec. VI, we will show our results for the 
two-dimensional Hubbard model. Comparisons between 
the DCA and the results of finite-sized simulations will 
be made in order to outline the complementarity of the 
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FIG. 1. A single step illustration of coarse-graining in the 
DMFA: all lattice propagators used to calculate the self energy 
are averaged over the points in the first Brillouin zone (top). 
This effectively maps the lattice problem to a single point in 
reciprocal space (bottom). Since the real space and reciprocal 
space of a single-site cluster are equivalent, this mapping takes 
the lattice problem to one of an impurity embedded within a 
host. 

MuUer-Hartmann ||] showed that this coarse-graining 
becomes exact in the limit of infinite-dimensions. For 
Hubbard-like models, the properties of the bare vertex 
are completely characterized by the Laue function A 
which expresses the momentum conservation at each ver- 
tex. In a conventional diagrammatic approach 



A(ki,k2,k3,k4 



J2 exp [ir ■ (ki -f k2 - kg - k4)] (2) 



ki+k2,k3-|-k4 



where ki and k.2 (ks and k4) are the momenta enter- 
ing (leaving) each vertex through its legs of G. However 
as the dimensionahty D oo Miiller-Hartmann showed 
that the Laue function reduces to M 



AD^oo(ki,k2,k3,k4) = 1 + 0{1/D) 



(3) 



The DMFA assumes the same Laue function, 
ADMFA(ki, k2, ka, k4) — 1, even in the context of fi- 
nite dimensions. Thus, the conservation of momentum 



at internal vertices is neglected. Therefore we may freely 
sum over the internal momentum labels of each Green 
function leg. This leads to a collapse of the momentum 
dependent contributions and only local terms remain. 

This argument may then be applied to the free energy 
functional. As discussed in many-body texts |]l9|, the 
additional free energy due to an interaction may be de- 
scribed by a sum over all closed connected graphs. These 
graphs may be further separated into compact and non- 
compact graphs. The compact graphs, which comprise 
the generating functional <&, consist of the sum over all 
single-particle irreducible graphs. The remaining graphs, 
comprise the non-compact part of the free energy. In the 
infinite-dimensional limit, $ consists of only local graphs, 
with non-local corrections of order 1/D. However, for the 
non-compact parts of the free energy, non-local correc- 
tions are of order one, so the local approximation applies 
only to <&. Thus, whereas irreducible quantities such as 
the self energy are momentum independent, the corre- 
sponding reducible quantities such as the lattice Green 
function are momentum dependent. 

The perturbative series for I] in the local approxima- 
tion is identical to that of the corresponding impurity 
model. However in order to avoid overcounting the lo- 
cal self-energy E(«a;„), it is necessary to exclude 'E{iujn), 
iujn = (2n -|- l)7rT is the Matsubara frequency, from the 
bare local propagator Q. Q{iLL)n)^^ — G{iujn)~^ + S(zw„) 
where G{iujn) is the full local Green function. Hence, 
in the local approximation, the Hubbard model has the 
same diagrammatic expansion as an Anderson impurity 
with a bare local propagator Q{iujn;'S) which is deter- 
mined self-consistently. 

An algorithm constructed from this approximation is 
the following: (i) An initial guess for E(za;„) is chosen 
(usually from perturbation theory), (ii) S(za;„) is used 
to calculate the corresponding local Green function 



G{iuj„ 



drj- 



(4) 



where p'^ is the non-interacting density of states, (iii) 
Starting from G{iujn) and S(iw„) used in the second step, 
the host Green function Q{iujn)~^ = G{iujn)~^ + S(iw„) 
is calculated which serves as bare Green function of the 
impurity model, (iv) starting with Q{iLUn), the local 
Green function GiitOn) is obtained using the Quantum 
Monte Carlo method (or another technique), (v) Using 
the QMC output for the cluster Green function G{iujn) 
and the host Green function Q^iun) from the third step, a 
new I](«w„) = Giii^n)^^ — G(iu!n)~^ is calculated, which 
is then used in step (ii) to reinitialize the process. Steps 
(ii) - (v) are repreated until convergence is reached. In 
step (iv) the QMC algorithm of Hirsch and Fye may 
be used to compute the local Green function G(r) or 
other physical quantities in imaginary time. Local dy- 
namical quantities are then calculated by analytically 
continuing the corresponding imaginary-time quantities 
using the Maximum-Entropy Method (MEM) ^ . 
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III. THE DYNAMICAL CLUSTER 
APPROXIMATION 

In this section, we will review the formalism which 
leads to the dynamical cluster approximation. Here, we 
first motivate the fundamental idea of the DCA which 
is coarse-graining, we then describe the mapping to an 
effective cluster problem and discuss the relationship be- 
tween the cluster and lattice at the one and two-particle 
level. 



A. Coarse-Graining 
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FIG. 2. A single step illustration of coarse-graining in the 
DCA: all lattice propagators used to calculate the self en- 
ergy are first averaged over the points within each cell in the 
Brillouin zone (top), mapping the lattice problem to a small 
cluster defined by the centers of the cells embedded within a 
host (bottom). 

Like the DMFA, the DCA may be intuitively motivated 
with a coarse-graining transformation. In the DMFA, the 
propagators used to calculate $ and its derivatives were 
coarse-grained over the entire Brillouin zone, leading to 
local (momentum independent) irreducible quantities. In 



the DCA, we wish to relax this condition, and systemat- 
ically restore momentum conservation and non-local cor- 
rections. Thus, in the DCA, the reciprocal space of the 
lattice (Fig. H) which contains N points is divided into 
Nc cells of identical linear size Ak. The coarse-graining 
transformation is set by averaging the Green function 
within each cell. If A'c = 1 the original lattice prob- 
lem is mapped to an impurity problem (DMFA). If Nc 
is larger than one, then non-local corrections of length 
« n/Ak to the DMFA are introduced. Provided that 
the propagators are sufficiently weakly momentum de- 
pendent, this is a good approximation. If Nc is chosen 
to be small, the cluster problem can be solved using con- 
ventional techniques such as QMC, NCA or FLEX. This 
averaging process also establishes a relationship between 
the systems of size N and Nc . A simp le and unique choice 
which will be discussed in Sec. [II B is to equate the ir- 
reducible quantities (self energy, irreducible vertices) of 
the cluster to those in the lattice. 



B. A diagrammatic derivation 

This coarse graining procedure and the relationship of 
the DCA to the DMFA is illustrated by a microscopic 
diagrammatic derivation of the DCA. 




FIG. 3. Coarse-graining cells for Nc = 8 (differentiated by 
alternating fill patterns) that partition the first Brillouin Zone 
(dashed line). Each cell is centered on a cluster momentum 
K (filled circles). To construct the DCA cluster, we map a 
generic momentum in the zone such as k to the nearest cluster 
point K = M(k) so that k = k — K remains in the cell around 
K. 

The DCA systematically restores the momentum con- 
servation at internal vertices relinquished by the DMFA. 
The Brillouin-zone is divided into Nc = cells of size 
Afc ^ 2tt/L (c.f. Fig. I for iV^ = 8). Each cell is rep- 
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resented by a cluster momentum K in the center of the 
cell. We require that momentum conservation is (par- 
tially) observed for momentum transfers between cells, 
i.e., for momentum transfers larger than A/c, but ne- 
glected for momentum transfers within a cell, i.e., less 
than Afc. This requirement can be established by using 
the Laue function B 



ADCA(ki,k2,k3,k4 



M(ki)+M(k2),M(k3)+M(k4 



connected compact diagrams, such as the one shown in 
Fig. ^. The corresponding DCA estimate for the free 
energy is 

Fdca = -ksT ($e - tr [E^G,] - tr In [-G,]) (7) 

where $c is the cluster generating functional. The trace 
indicates summation over frequency, momentum and 
spin. Foe A is stationary with respect to G^, 



(5) 



where M(k) is a function which maps k onto the momen- 
tum label K of the cell containing k (see. Fig. ||). This 
choice for the Laue function systematically interpolates 
between the exact result, Eq. 0, which it recovers when 
Nc ^ N and the DMFA result, Eq. ^, which it recovers 
when Nc — I. With this choice of the Laue function the 
momenta of each internal leg may be freely summed over 
the cell. 

This is illustrated for the second-order term in the gen- 
erating functional in Fig. ^. Each internal leg G(k) in a 
diagram is replaced by the coarse-grained Green function 
G(M(k)), defined by 



N 



(6) 



where N is the number of points of the lattice, Nc is the 
number of cluster K points, and the k summation runs 
over the momenta of the cell about the cluster momentum 
K (see. Fig. |^). The diagrammatic sequences for the 
generating functional and its derivatives are unchanged; 
however, the complexity of the problem is greatly reduced 
since Nc <^ N. 




K+Q 



^DCA 



N 



G(K) 




TG(KH-k): 

^ K'+Q~ 

FIG. 4. A second-order term in the generating functional 
of the Hubbard model. Here the undulating line represents 
the interaction U, and on the LHS (RHS) the solid line the 
lattice (coarse-grained) single-particle Green functions. When 
the DCA Laue function is used to describe momentum conser- 
vation at the internal vertices, the momenta collapse onto the 
cluster momenta and each lattice Green function is replaced 
by the coarse-grained result. 

As with the DMFA, the coarse-graining approxima- 
tion will be applied to only the compact part of the free 
energy, $, and its derivatives. This is justified by the 
fact that there is no need to coarse-grain the remaining 
terms in the free energy. Formally, we have justified this 
elsewhere by exploring the Afc-dependence of the com- 
pact and non-compact parts of the free energy pO| . The 
generating functional is the sum over all of the closed 



-1 SFpcA 



Ee.(M(k)) - E,(k) = 0, 



(8) 



which means that S(k) = Eco.(M(k)) is the proper ap- 
proximation for the lattice self energy corresponding to 
The corresponding lattice single-particle propagator 
is then given by 



G(k,z) = 



1 



z-ek-e-S,(M(k),z) 



(9) 



A variety of techniques may be used to sum the cluster 
diagrams in order to calculate Ec and the vertex functions 
Tc- In the past, we have used QMC [^, the non-crossing 
approximation |^ or the Fluctuation-Exchange approxi- 
mation. Here, we will use the QMC technique which we 
will detail in Sec. IV. Since QMC is systematically ex- 



act; i.e. it effectively sums all diagrams to all orders, care 
must be taken when defining the initial Green function 
(the solid lines in Fig. ^ to avoid overcounting diagrams 
on the cluster. For example, to fourth order and higher 
in perturbation theory for the self energy, non-trivial self 
energy corrections enter in the diagrammatic expansion 
for the cluster self energy of the Hubbard model. To avoid 
overcounting these contributions, we must first subtract 
off the self energy corrections on the cluster from the 
Green function line used to calculate Ec and its func- 
tional derivatives. This cluster-excluded Green function 
is given by 



1 



1 



g{K,z) G(K,z) 



Sc(K,z) 



(10) 



which is the coarse-grained Green function with corre- 
lations on the cluster excluded. Since Ec(K,z) is not 
known a priori, it must be determined self-consistently, 
starting from an initial guess, usually from perturbation 
theory. This guess is used to calculate G from Eq. ||. 
g(K, z) is then calculated with Eq. and it is used to 
initialize the QMC calculation. The QMC estimate for 
the cluster self energy is then used to calculate a new 
estimate for G(K) using Eq. ^. The corresponding Q{K.) 
is used to reinitialize the procedure which continues un- 
til Gc = G and the self energy converges to the desired 
accuracy. 

One of the difhculties encountered in earlier attempts 
to include non-local corrections to the DMFA was that 
these methods were not causal [^l]j2^. The spectral 
weight was not conserved and the imaginary parts of the 
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one-particle retarded Green functions and self-energies 
were not negative definite as required by causality. The 
DCA algorithm presented in this subsection does not 
present these problems. This algorithm is fully causal 
as shown by Hettler et al. |^ . They analyze the different 
steps of the self-consistent loop and found that none of 
them breaks the causality of the Green functions. Start- 
ing from the QMC block, one can see that if the input 
Q is causal, since the QMC algorithm is essentially ex- 
act, the output Gc will also be causal. Then the cor- 
responding I]c(K,zcij„) is causal. This in turn ensures 
that the coarse-grained Green function G(K,ia;„) also 
fulfills causality. The only non-trivial operation which 
may break causality is the calculation of tJ(K, iuJn)- Het- 
tler et al. used a geometric proof to show that even this 
part of the loop respects causality. 

In the remainder of this section, we will give further 
details about the DCA formalism, and discuss the re- 
lationship between the cluster and the lattice problems. 
Below, we will discuss the steps necessary to choose the 
coarse-graining cells and ensure that symmetries of the 
lattice are preserved. 



C. Selecting the Coarse-Graining cells 
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As we will see in Sec. IV the solution of the cluster 



problem using the quantum Monte Carlo method though 
is a great simplification over the original lattice prob- 
lem is still a formidable task. The reason is that the 
self-consistent nature of the cluster problem forces us to 
adopt the Hirsch-Fye algorithm. While this algorithm 
is very efficient for few impurity problems, it becomes 
slow even for a cluster of a modest size. Therefore, in 
order to study the size dependence of physical quanti- 
ties we adopt various cluster tilings of the lattice instead 
of confining ourselves to only the usual square tilings 
A^c^ 4, 16,36,64.... 
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FIG. 5. Different tile sizes and orientations in real space. 
The tiling principle translation vectors, ai and a2, form two 
sides of each tiling square (illustrated for the Nc = 20 tiling). 
For square tile geometries, a2x — —aiy and a2y — aix. 

When selecting the coarse-graining cells, it is impor- 
tant to preserve the point group symmetries of the lat- 
tice. For example, in this study, we will choose a simple 
square lattice. Both it and its reciprocal lattice share 
6*41, symmetry with eight point group operations. We 
must choose a set of coarse-graining cells which preserve 
the lattice symmetry. This may be done by tiling the 
real lattice with squares, and using the K points that 
correspond to the reciprocal space of the tiling centers. 
We also will only consider tilings which contain an even 
number of sites, to avoid frustrating the magnetic correla- 
tions on the cluster. Square tilings with an even number 
of sites include = 4,8,10,16,18,20,26,32,34,36,.... 
The first few are illustrated in Fig. The relation be- 
tween the principle lattice vectors of the lattice centers, 
ai and a2 , and the reciprocal lattice takes the usual form 
= 27raj/(ai x a2), with K„„ = ngi -I- mg2 for inte- 
ger n and m. For tilings with either aix = aiy (corre- 
sponding to Nc — 1, 8, 18, 32...) or one of aix or aiy zero 
(corresponding to Nc — 1, 4, 16, 36...), the principle recip- 
rocal lattice vectors of the coarse-grained system either 
point along the same directions as the principle recipro- 
cal lattice vectors of the real system or are rotated from 
them by 7r/4. As a result, equivalent momenta k are al- 
ways mapped to equivalent coarse-grained momenta K. 
An example for Nc = 8 is shown in Fig. ^. However, 
for Nc — 10,20,26,34..., the principle reciprocal lattice 
vectors of the coarse-grained system do not point along 
a high symmetry direction of the real lattice. Since all 
points within a coarse-grained cell are mapped to its cen- 
ter K, this means that these coarse-graining choices vio- 
late the point group symmetry of the real system. This 
is illustrated for A^c = 10 in Fig. ^, where the two open 



6 



dots resting at equivalent points in the real lattice, fall in 
inequivalent coarse-graining cells and so are mapped to 
inequivalent K points. Thus the tilings corresponding to 
Nc = 10, 20, 26, 34... violate the point-group symmetry of 
the real lattice and will be avoided in this study. 



Nc=10 





FIG. 6. The coarse graining cells for Nc — 8 and 10 each 
centered on a coarse-grained momenta K represented as black 
filled dots. For Nc = 8 equivalent momenta k are always 
mapped to equivalent coarse-grained momenta K. However, 
this is not true for Nc = 10 where, for example, the two 
equivalent momenta shown by open dots are mapped to in- 
equivalent coarse-grained momenta. 

One should note that the coarse-graining scheme also 
depends strongly on dimensionality. For example, in one 
dimension, any cell with an even number of sites will 
preserve the lattice symmetry. 



D. One-particle Green functions 

In the DMFA, after convergence, the local Green func- 
tion of the lattice is identical to that of the impurity 
model. Though in the DCA, the coarse-grained Green 
function G(K, iujn) is equal to the cluster Green function 
Gc(K,za;„), this quantity is not, however, used as an ap- 
proximation to the true lattice Green function G(K, iuJn). 
The correct procedure to calculate the lattice physical 
quantities within the DCA is to approximate the lat- 
tice irreducible quantities with those of the cluster. The 
lattice reducible quantities are then deduced from the 
irreducible. This procedure was justified formally in 



Sec. Ill B. To obtain a physical understanding, one must 
first understand why reducible and irreducible quanti- 
ties must be treated differently. Consider a quasiparticle 
propagating through the system. The screening cloud is 
described by the single-particle self-energy S(k, lu) which 
itself may be considered a functional of the interaction 
strength U and the single-particle propagator G(k, w), 
S = E[J7, G]. The different screening processes are de- 
scribed perturbatively by a sum of self-energy diagrams. 
If the size of the screening cloud rg is short, the propaga- 
tors which describe these processes need only be accurate 
for distances < rg. From the Fourier uncertainty prin- 
ciple, we know that the propagators at short distances 
may be accurately described by a coarse sampling of the 
reciprocal space, with sampling rate Ak = Tr/rg. Hence, 
in this case, S[J7, G] may be quite well approximated by 



S[f/,G]. 

On the other hand, the phase accumulated as the par- 
ticle propagates through the system is described by the 
Fourier transform of the single-particle Green function. 
Since this accumulated phase is crucial in the descrip- 
tion of the quantum dynamics it is important that G(r) 
remains accurate at long distances, so it should not be 
coarse-grained as described above. However it may be 
constructed from the approximate self-energy. Hence, 
the approximate lattice Green function is given by Eq. ||. 
Thus, as in the DMFA, the lattice Green function is gen- 
erally more strongly momentum dependent than the cor- 
responding self energy. 

In the case of the 2D Hubbard model, non-local corre- 
lations are the most important in the parameter regime 
close to the quantum critical point at half filling. Away 
from this parameter regime rg is thus expected to be 
short. Here, the above construction scheme for the ap- 
proximate lattice Green function is likely to yield accu- 
rate results even for clusters of modest size. However, as 
the quantum critical point is approached, longer range 
correlations are important. As a consequence one will 
need to evaluate S[?7, G] on larger clusters. 



E. Two-Particle Green Functions 

A similar procedure is used to construct the two- 
particle quantities needed to determine the phase dia- 
gram or the nature of the dominant fluctuations that can 
eventually destroy the quasi-particle. This procedure is 
a generalization of the method of calculating response 
functions in the DMFA l23|j2l. In the DCA, the intro- 



duction of the momentum dependence in the self-energy 
will allow one to detect some precursor effects which are 
absent in the DMFA; but for the actual determination 
of the nature of the instability, one needs to compute 
the response functions. These susceptibilities are ther- 
modynamically defined as second derivatives of the free 
energy with respect to external fields. ^c{Gr) and Eco-, 
and hence Foe A depend on these fields only through Go- 
and G°. Following Baym it is easy to verify that, 
the approximation 



- ca,a' 



(11) 



yields the same estimate that would be obtained from 
the second derivative of FucA with respect to the ap- 
plied field. For example, the first derivative of the free 
energy with respect to a spatially homogeneous external 
magnetic field h is the magnetization, 

m = tr[crG^]. (12) 

The susceptibility is given by the second derivative. 



dm 
'dh 



= tr 



dG„ 



dh 



(13) 
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We substitute 
derivative, 



dm 
'dh 



= tr 



dG„ 



dh 



G 



^ tr 



0-1 



GMl 



and evaluate the 



dG^> dh 



(14) 

If we identify Xo-,<t' = g ^g^' , and x° = G^, collect all 
of the terms within both traces, and sum over the cell 
momenta k, we obtain the two-particle Dyson's equation 



mo 







(15) 



We see that again it is the irreducible quantity, i.e. the 
vertex function, for which cluster and lattice correspond. 



1. particle-hole 



By coarse-graining the Bethe-Salpeter equation, we have 
greatly reduced its complexity; each of the matrices 
above is sufficiently small that they may be easily ma- 
nipulated using standard techniques. 

In contrast with the single-particle case where the 
coarse-grained quantities are identical to those of the 
cluster, Xca,a'{'i^K,K') is not equal to Xcr,<7'{q, K, K'). 
This is because the self-consistency is made only at 
the single-particle level. Unlike the single particle case 
where both T,{K) and G{K) are directly calculated, nei- 
ther T„ ,ji{q,K,K') nor the coarse-grained susceptibility 
X<t,ct' {q, K, K') are calculated during the self-consistency. 
Instead, the coarse-grained non-interacting susceptibility 
X° is calculated in a separate program after 

the DCA converges using the following relation 



N 



^ G,(K + k, ic^„)G,(K + k + q. 



(18) 



In this subsection we will provide more details about 
the relationship between the lattice and cluster two- 
particle Green functions and describe how a particle- hole 
susceptibility may be calculated efficiently. As a spe- 
cific example, we will describe the calculation of the two- 
particle Green function 

j-P rf3 t-P 

Xcr,a'{q,k,k') ^ / / / / dTldT2dT3dT4 

Jo Jo Jo Jo 

^ gj(([i;„+i^)ri-cj„r2+cj„/r3-(w„/+;y)r4) 

^ (rrCk+q^('ri)Ck^('^2)ci,,^,(r3)Ck,+q<,,(r4)) 

where we adopt the conventional notation [|l9| k = 
(k, iwn), k' = (k, w^), q = (q, i^n) and Tr is the time 
ordering operator. 

Xa-,(T'{q,k,k') and Tcy.a-'iq,k,k') are related to each 
other through the Bethe-Salpeter equation: 

Xa,a' {q, k, k') = X°a,a' k, k') + xl,a" (q, k, k") 

X T,„^,,„{qX,k"')xa"'.a'{q,k"',k') (16) 

where T„^fji(q, k, k') is the two-particle irreducible vertex 
which is the analogue of the self-energy, x° {q, k, k") is 
the non-interacting susceptibility constructed from a pair 
of fully-dressed single-particle Green functions. As usual, 
a summation is to be made for repeated indices. 

We now make the DCA substitution Fo-.o-' (q, k, k') 
re<^,„/ (q,M(k),M(k')) in Eq. |l6| (where ' frequency la- 
bels have been suppressed) . Note that only the bare and 
dressed two-particle Green functions x depend upon the 
momenta k within a cell. Since x and x*^ in the prod- 
uct on the RHS of Eq. |l^ share no common momentum 
labels, we may freely sum over the momenta k within a 
cell, yielding 

Xa.a' (g, K') = xl,a' {q. K,K')+ xl^,„ (g, K, K") 
X T,,„^„,„{q,K\K"')x.'",.'{q.K"',K'). (17) 



The corresponding cluster susceptibili ty is calculated in 
the QMC process, as discussed in Sec. IVD and the ver- 
tex function is extracted by inverting the cluster two- 
particle Bethe-Salpeter equation 

Xc,,,, (g, K, K') = xcl.a' {q, K, K') + Xc°,,. (g, K, K") 

xr,,„^,,„ (g, K\ K"')xcan,,„, {q, K'" , K') . (19) 



If we combine Eqs. 19 and |l^, then the coarse-grained 
susceptibility may be obtained after elimination of 
V{q,K,K') between the two equations. It reads 



(20) 



where, for example, x is the matrix formed from 
Xa.a' {q, K, K') for fixed q. The charge (c/i) and spin [sp) 
susceptibilities Xch.sp{q,T) are deduced from x 



Xch.sp{q-,T) 



X! ^'y'y'Xa,a'{q,K,K') 



(21) 



if A'' 



where Ao-o-' — 1 for the charge channel and A^ct' = cnr' 
for the spin channel. 



2. particle-particle 

The calculation of susceptibilities in the particle- 
particle channel is essentially identical to the above. The 
exception to this rule occurs when we calculate suscepti- 
bilities for transitions to states of lower symmetry than 
the lattice symmetry. For example, in order to obtain 
the pair function of the desired symmetry (s,p,d), the 
two-particle Green function must be multiplied by the 
corresponding form factors (?(k) and (?(k'). In the study 
of the Hubbard model below, we will be particularly in- 
terested in g(k) — 1 {s wave), g(k) = cos{kx) + cos{ky) 
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(extended s wave) and .g(k) = cos{kx) — cos{ky) (dx^-y^ 
wave). These symmetries have been evoked as possible 
candidates for the superconducting ground state. 
These factors modify the Bethe-Salpeter equations 

5(k)x(g, fc, fc')5(k') = 5(k)x°(9, k, k')g{\^') (22) 
+ .9(k)x°(g, fc, k") X r(g, fc", k'") X x(9, fc'", A;')g(k') 



where 
Xi<l,k,k') 



X e 



/3 ./3 ./3 ./3 

/ / / dTldT2dT^dTi 

Jo Jo Jo 



(23) 



F. Local quantities 

We will also need to evaluate a number of local quan- 
tities on the lattice. They include the magnetic moment, 
the local magnetic susceptibility, the local Green func- 
tion, etc. The local cluster quantities are identical to the 
local lattice ones. This may be seen for example on the 
one-particle Green function. The coarse-grained Green 
function is related to the lattice Green function as fol- 
lows 



^(^' ^) = ^ E E e'^K-('--')e^*^-(^+'-')G(X + r', lo). (31) 



K,k- 



X (7V4+q^(^l)cLk-o-(^2)C-k'-<T(T3)Ck'+q^(T4)) 

On the LHS, we have dropped the spin indices since we 
will consider only opposite-spin pairing. Eq. ^ cannot be 
easily solved if it is coarse-grained, since this will partially 
convolve x(q, fc, k') with two factors of g on the LHS and 
one factor on the RHS. Hence for the pairing susceptibil- 
ities, or for any situation where non-trivial form factors 
must be used, we use the equivalent equation involving 
the reducible vertex T2 (instead of the irreducible vertex 

r) 

.g(k)x(g, fc, fc').9(k') - .9(k)x"(g, k, k')g{V!) 

+ 5(k)x^g,fc,fc") 

xT2{qXX')x\lX'.k')g{V!), (24) 

where 

T2{q,k,k')^T{q,k,k') (25) 
+ x"(9,fc,A:")r(g,fc",A"')x°(9,fc"',fc') + --- 



We define 



Hg,g {q,k,k')^ .9(k)x(g, fc, fc').9(k') (26) 
H°_g((Z, fc, k') = g{\^)x°{q. k, k')g{k') (27) 
U"g{q,k,k')=g{-k)x"{q,k,k'). (28) 

The remaining steps of the calculation are similar to the 
particle-hole case. We invert the cluster particle-particle 
Bethe-Salpeter equation with g = I for the cluster, in 
order to extract Tc- We then coarse-grain Eq. |2^, and use 

Tc to calculate the coarse-grained T2 = Fc (l — X^^c) 
We then coarse-grain Eq. |2^, and use the coarse-grained 
T2 to calculate the coarse-grained flg,g 



It is easy to see from this relation that G{0, uj) = G(0, lu) 



IV. THE QUANTUM MONTE CARLO 
ALGORITHM 

In this section we will derive a generalization of 
the Hirsch-Fye Anderson impurity algorithm suitable 
to simulate a Hubbard cluster embedded in a self- 
consistently determined host. We will then discuss the 
differences between this algorithm and the more famil- 
iar Blanckenbecler-Sugar-Scalapino (BSS) algorithm p^ ] 
used to simulate finite-sized systems. Finally, we will 
discuss how different quantities mentioned above may be 
measured efficiently and how the code can be optimized. 



A. Formalism 

The Hirsch-Fye algorithm is an action-based technique. 
Therefore, knowledge of the underlying Hamiltonian is 
not required provided that we know the Green function 
for the non- interacting cluster coupled to the host, and 
the interacting part of the action or Hamiltonian. The 
interacting part is unchanged by the coarse-graining since 
it is purely local. It may be written in the real space as 
follows, 



i/7 = C/^(ni,T-l/2)(ni,x-l/2), 



(32) 



i=l 



Ilg^g{q,K,K')^IllJq,K,K' 



(29) 



U%, K, K")T2{q, K", K"')UUq, K'", K') 



The pairing susceptibility of a desired symmetry is given 

by 



^.(9'^) = ^^En..('z,if,^') 



and the bare cluster Green function is ^?(K, zw„). 

Given this information, the most direct way to derive 
this algorithm is to express the partition function as path 
integrals over Grassmann variables. The first step is to 
disentangle the interacting, iJ/, and non-interacting, Hq, 
parts of the Hamiltonian using a Trotter-Suzuki decom- 
position for the partition function. We divide the interval 
[0,/3] into Ni sufficiently small subintervals At — /3/Ni 
such that Ar^ [Ho, Hj] may be neglected. This leads to 



(30) 



N, 



N, 



K.K' 



-AtH 



TrY[, 



-AtHo^-AtHi 



(33) 



1=1 
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The interacting part of the Hamiltonian may be further 
decoupled by mapping it to an auxihary Ising field via a 
discrete Hirsch-Hubbard-Stratonovich (HHS) ||2^ trans- 
formation, 



ArHj ^ g-AT(7^.(ni-r-l/2)(nii-l/2) 



e ' = e 

^ lg-AT[//4 

2 



HE' 

i Si=±l 



Si(niT 



(34) 



where cosh(a) = e^'^'^^^. 

We now introduce coherent states of the operators on 
the cluster and in the host as the basis states and ex- 
press the partition function as path integrals over the 
corresponding Grassmann variables 7i,i,CT and (/'k.i.o- de- 
fined over each TV/ time slices r/ = /Ar of the interval 
[0,/3] ||2|]. After substituting the Grassmann variables 
one obtains the following approximation for the partition 
function which becomes exact as At — > 0, 



Z 



-So[7.0]g-S/[7] 



(35) 



where symbols denote the measures of path integra- 
tion over the corresponding Grassmann fields and S'(o)/ is 
the (non-)interacting part of the action. The interacting 
part of the action, becomes 



A'c Ni 

•^/W = - E E E "7M,^crsi,i7i. 

i=l 1 = 1 a 



l-l,<y 



(36) 



The non-interacting parts are 

/ 0k,i. 



5o[7,^] = At^ 



Ok.i-l., 



i.l.cr 



a7 



(37) 



since, as we show below, we only require knowledge of the 
ratio of the partition functions for two different config- 
urations of the HHS fields. Other fixed prefactors (de- 
pending upon U, (3 ■ • ■) have also been disregarded in 
Eq. pa. Sc is the cluster action. It takes the form 



Sc[l] = 



7^.,g-i(i,Z;i',Z')7i'.P,^ + 5/[7] (39) 



where CJ(i, /; i', Z') is the cluster excluded (i.e. non- 
interacting on the cluster) Green function, defined pre- 
viously. Now we will integrate out the remaining cluster 
Grassmann variables. The partition function then be- 
comes 



(40) 



where again factors which are fixed during the QMC pro- 
cess have been ignored. {Gca-sii) 
Green function matrix with elements 



is the inverse cluster 



-1 



5i,j<^/',i-iao-Si(r;) + Q-^ju 



(41) 



If we re-exponentiate the first term in the RHS of the 
above formula by defining Va-(i, = otSi^ia, we can write 
Eq. P in a simple matrix notation as 



T e^' - 1 



(42) 



where T is Sij6i-i^i'. The matrix product Gc^ e " de- 
pends upon the HHS fields only along its diagonal ele- 
ments. As can be seen from Eq. 37, each diagonal ele- 
ment of the matrices and hence Gc^^ is 1. Therefore, 
the inverse Green functions for two different field config- 
urations, {si;} and are related by 



Gl-'e-^'^ 



G.-'e-^' 
Or, after multiplying by " , 



and collecting terms 



At H^ciusteMl,i,aAKi-iynli,ani,i-i,.) G,'/ - G,-' = (G,-' - l)e-^"(e 



V' 



(43) 



(44) 



k,i,/,cr 

where Hhost, H^iuster Hamiltonian for the host, 

and the non-interacting degrees of freedom on the cluster 
including the coupling to the host, respectively. The de- 
tailed form of both Hhost and H^^^^^^^ are unknown, due 
to the self-consistent renormalization of the host. How- 
ever, both are purely bilinear, and may be integrated out 
of the action without further approximation. 

We will first integrate out the host degrees of freedom. 
The partition function becomes 



Multiplying from the left by Gc and from the right by 
Gc', we find 



Gc n 



or 



f ^ 



Gc<. + (Gc, -l)(e^^-^'-l)G^ 
l + (l-G,.)(e^^-^"-l). 



(45) 
(46) 



Z cx 



]Jdet(gk,^) ^ 



(38) 



where gk.cr is the Green function of the host. It remains 
fixed during the QMC process, and it may be disregarded 



B. The QMC algorithm 

We will now proceed to derive the Monte Carlo al- 
gorithm. The QMC algorithm involves changes in 
the Hubbard-Stratonovich field configuration {si,;} — > 
{s'i,;}, and accepts these changes with the transition 
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probability Ps—>s' ■ Thus, to define the algorithm, we need 
Ps^s' and a relation between the cluster Green functions 
G and G" for the two different auxiliary field configura- 
tions. To simplify the notation, we introduce a combined 
space-time index i — (i, Z), and will consider only local 
changes in the fields s,„ ~* s' „i = —s„i. As can be in- 
ferred from Eq. the probability of a configuration {si] 



is Ps oc det(G'c^j^})det(Gc^{^^.j)^ 
tailed balance requires Ps'Pg'^s = PsPs 



on the other hand de- 



for all 



We may satisfy this requirement either by defining the 
transition probability Ps'^s = R/{^ + R), where 



R = 



Ps' ^ det(GcT)det(Gex) 



(47) 



is the relative weight of two configurations, or by letting 
Ps'->s = minimum (i?, 1) (the first choice is called the 
"heat bath" algorithm, and the second the "Metropolis" 
algorithm). If the difference between two configuration 
is due to a flip of a single Hubbard Stratonovich field at 
the mth location in the cluster space-time , then from 
Eq. M 



R = l[[l + {1-Gc 



J(e- 



-aCT(s„ 



1)]" 



(48) 



For either the Metropolis or the heat bath algorithm, if 
the change is accepted, then we must update the Green 
function accordingly. The relationship between G and G' 
is defined by Eq. |4^ 



1 + (1 - G,„)(e-"-(^'"--'-) _ 1) 



G, 



camj • 



(49) 



The QMC procedure is initialized by setting Gcaij — 
Qij where Qij is the (cluster) Fourier transform of 5(K) 
(Eq. |l^, and choosing the correspondin g fie ld configu- 
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to create a 



ration with all Si = 0. Then we use Eq. 
Green function corresponding to a meaningful field con- 
figuration (i.e. Si = ±1, for each i — or the {si} 
from a previous run or iteration). We proceed by se- 
quentially stepping through the space-time of the clus- 
ter, proposing local changes Si — > — s^. We accept the 
change if Pg'^s is greater than a random number between 
zero and one and update the Green function according to 
Eq. (^). After twenty to one hundred warm-up sweeps 
through the space-time lattice of the cluster, the sys- 
tem generally comes into equilibrium and we begin to 
make measurements. A few lattice updates are used be- 
tween each measurement step to reduce the correlations 
between measurements. This improves the efficiency of 
the algorithm, since as we will see below, the measure- 
ments are numerically expensive. After many iterations 
of lattice updates, numerical round-off error begins to ac- 
cumulate in the Green function update, Eq. ^ To com- 
pensate for this round-off error, the Green functions must 
be refreshed by again setting Gc^ij — Qij , and then using 
Eq. ^ to recalculate the Green function corresponding 
to the present field configuration. 



C. Differences with the BSS Algorithm 

The Hirsch-Fye (HF) algorithm differs in several ways 
from the more familiar Blanckenbecler-Sugar-Scalapino 
(BSS) algorithm |26| used to simulate finite-sized sys- 
tems. 

The BSS algorithm is more efficient. HF simulations 
can be computationally quite expensive since the mem- 
ory and the CPU time required by this algorithm scale as 
(iVciV;)^ and {NcNi)^, where Nc and Ni are respectively 
the number of cluster sites and the number of time slices. 
The BSS algorithm scales as NiN^ for the memory and 
NiN^ for the CPU time. In order to study a meaningful 
set of cluster sizes using the Hirsch-Fye algorithm, it is 
necessary to use massively parallel computers. The max- 
imum size we studied is = 64 for the two-dimensional 
Hubbard model. This maximum size is indeed smaller 
than what can be reached with the BSS algorithm ap- 
plied for finite system simulations (FSS). But, one should 
bear in mind that, in the DCA, the system is in the ther- 
modynamic limit, it is the range of spatial correlations 
which is restricted to the cluster size. Cluster size ef- 
fects are of different nature than that occurring in FSS. 
Therefore, the DCA as discussed in previous studies, can 
provide information which cannot be obtained from the 
FSS. 

The Hirsch-Fye algorithm is action-based, whereas the 
BSS algorithm is Hamiltonian based. Therefore, the BSS 
algorithm cannot be employed to solve the DCA cluster 
problem, since the cluster problem has no Hamiltonian 
formulation with known parameters, and its action is 
highly non-local in time. The BSS algorithm requires 
that the action be local in time. The cluster action, 
Eq. |3^, is long-ranged in time due to the term involving 
Q. Thus, the Hirsch-Fye algorithm is the most appropri- 
ate QMC algorithm to solve the DCA embedded cluster 
problem. 

In addition to the differences mentioned above (de- 
tailed knowledge of the Hamiltonian is not needed for 
the HF algorithm so long as we have an initial Green 
function), there are other advantages to the HF algo- 
rithm. Whereas in the BSS algorithm, all degrees of 
freedom must appear explicitly, in the HF algorithm, any 
non-interacting degrees of freedom may be integrated out 
without loss of any information. At the end of the cal- 
culation, the irreducible diagrams on the interacting or- 
bitals may be used to calculate any relevant quantity. 
Therefore, the HF algorithm may be used, for example, 
to simulate the periodic Anderson model (with only the 
f-orbital correlated) with the same computational cost 
required to simulate a single-band model. One may also 
incorporate a (dynamical) mean field coupling to an en- 
vironment, or between an infinite set of coupled Hubbard 
planes at no additional computational cost. In these 
cases, the information about the mean field coupling to 
the environment or the other planes is reflected in Q. 

For clusters, the Hirsch-Fye algorithm is very stable at 
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low temperatures. In particular, the matrices which 
are generated in the algorithm are well-conditioned, the 
costly stabilization steps required at low temperatures 
po| , pi| for the more popular BSS algorithm are avoided. 

Finally, the Hirsch-Fye algorithm is easily adapted to 
making measurements which are non-local in time, such 
as those required to calculate the irreducible vertex func- 
tions. This will be discussed in the next section. It is very 
difficult to measure quantities which are non-local in time 
with the BSS algorithm. In fact, such measurements re- 
quire significantly more CPU time than is required to av- 
erage over the HHS field configurations, since the CPU 
time required by these measurements scales like {NcNiY 
for both the BSS and Hirsch-Fye algorithms. Thus, when 
these measurements which are non-local in time are re- 
quired, both algorithms scale like [NcNiY . 



D. Making and Conditioning Measurements 

In the QMC technique, all the physical quantities are 
expressed in terms of Green functions. Standard dia- 
grammatic techniques are applied to evaluate these quan- 
tities. In doing so one must remember that the Hubbard- 
Stratonovich transformation reduces the problem to one 
of free electrons moving in a time-dependent field. Thus 
for each field configuration, any diagram may be formed 
by summing all allowed Wick's contractions. The full 
quantity is recovered by averaging this over all field con- 
figurations. Connected as well as disconnected configura- 
tions must be used during the evaluation. It is important 
to average over all equivalent time and space differences 
and all the symmetries of the Hamiltonian in order to 
produce the lowest variance measurement. 

One difficulty encountered with the DCA algorithm 
is that a reliable transform from imaginary-time quan- 
tities, in the QMC part, to Matsubara frequencies, for 
the coarse-graining part is needed. A careful treatment 
of the frequency summation or the imaginary-time in- 
tegration is crucial in order to ensure the accuracy and 
the stability of the algorithm and to maintain the correct 
high-frequency behavior of the Green functions. We need 
to evaluate the following integral 



Gc (K,ia;„) 



dre*'^"^Gc(K,T) 



(50) 



But from the QMC, we know the function Gc(K, r) only 
at a discrete subset of the interval [0, /?]. As it may 
be readily seen by discretizing the above equation, the 
estimation of Gc (K,ja;„) becomes inaccurate at high- 
frequencies. This is formalized by Nyquist's theorem 
which tells us that above the frequency — unpre- 
dictable results are produced by conventional quadrature 
techniques. For example, a rectangular approximation to 
the integral in Eq. ^yields a G(K,ia;„) that is periodic 
in ujn- This presents a difficulty since the causality re- 
quires that 



lim G(K, iLOn) ~ -. — . 

;„-+oo iu) 



(51) 



A straightforward way to cure this problem may be to 
increase the size of the set of r-points where the Green 
function is evaluated. But, this renders the QMC simu- 
lation rapidly intractable as seen in the previous section. 
A more economic way to avoid the problem is to use the 
high frequency information provided by an approximate 
method that is asymptotically exact. 

Second-order perturbation theory is enough to obtain 
the correct asymptotic behavior, Eq. To use this 
high frequency information, we compute the Matsubara- 
frequency Green function from the imaginary-time QMC 
Green function as follows 1321 



Gc(K, zcj„) =Gcpt(K,iLj, 



f dTe*'^^(G,(K,r)-G,pt(K,r)) 
Jo 



(52) 



The integral is computed by first splining the difference 
Gc(K, r) — Gcpt(K, r) using an Akima sphne Q, and 
then integrating the spline (a technique often called over- 
sampling) . 

As another example, consider the local magnetic sus- 
ceptibility (used to calculate the screened local moment) 

^-^iijy^ wciiM4i(o)^Ti(o)) 



xG, 



r(i,T';i,T + T')){,,.j 



(53) 



where the {su] subscript indicates that the Monte Carlo 
average over the Hirsch-Hubbard-Stratonovich fields is 
still to be performed, and in the last step in Eq. ^ we 
form all allowed Wick's contractions and average over 
all equivalent time and spatial differences to reduce the 
variance of this estimator. This measurement is best ac- 
complished by splitting it in two parts. First, we measure 
X(r) 



xG,_,(i,r';i,r + r')){,,,} (54) 



by approximating the integral as a sum using a rectan- 
gular approximation. For r > 



xin) 



2NiN, 



xGc_,(i,r;i,znrf(/ + ?'))){si,} 



(55) 
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where ind(l) is the smaller nonnegative value of either 
I or I — Ni. For r = the fact that we always store 
Gaih I'; i, I') — Gcr{i, + 0+; i, tv) requires us to modify 
the measurement 

x(-o).^E(Gc.(M';M') 

(7,1' ,i 

(G,_,(i,Z';i,n (56) 



Finally 



X(T)= rdTx(r)«^/(OATx(rO: 
Jo , 



(57) 



where the Simpson factor f{l) = 2At/3 (4Ar/3) for odd 
(even) Z is used to reduce the systematic error of the 
integral. 

As a final example, consider the cluster particle-particle 
Green function mat rix Xc (q, K, K') {K = (K,iti;„)) 
which is used in Sec. HIE 2 to calculate the lattice pair- 
field susceptibilities. The first step is to form the corre- 
sponding quantity in the cluster space-time 

Xc{Xi,X2,X3,Xi) = (TrC^iX,)ciiX2)cliX3)c\{X^) 



(58) 

Here Xi is in the space- (imaginary)time notation Xi ~ 
(Xi, Ti), where the points Xi are on the corresponding re- 
ciprocal cluster of K in real space, and (T,-...} denotes the 
time ordered averaging process. The two-particle Green 
functions are difficult to measure efficiently. For a partic- 
ular configuration of the auxiliary Hubbard-Stratonovich 
fields, the fermions are noninteracting, thus the expec- 
tation value indicated above may be evaluated in two 
steps. First, using Wick's theorem, its value is tabulated 
for each field configuration {si} and then transformed 
into the cluster Fourier space. Second, we Monte Carlo 
average over these configurations. After the first step, 
the expression for the above two-particle Green function 
in the cluster momentum-frequency space becomes 



Xc(Q, ii^n; K, iun; K', iw„ 



( e'''''''G,^{X,,X, 



4)e 



X\ ,X4 



E e^^Q~^'^^^G,,(X2,X3)e-W-^)^-^){,^j. (59) 

where K is the momentum-frequency point K — 
(Kjiojn)- The average over Hubbard-Stratonovich fields 
can be evaluated through the QMC sweeps along 
with the measurements of Gc^ and Gci- However, the 
sums (integrals) over r in Eq. ^ require special consid- 
eration. Since the Green functions change discontinu- 
ously when the two time arguments intersect, the best 
applicable integral approximation is the trapezoidal ap- 
proximation. Using this, we will run into Green func- 
tions G'c(X, r; X, t) with both time and space arguments 



the same. In the QMC algorithm, this is stored as 
Gc{'^,t~^;^,t) (i.e. it is assumed that the first time ar- 
gument is slightly greater than the second); however, 
if we replaced the equal time Green function to be 
the average {Gc(X,r+;X,r) Gc(X, t; X, r+)}/2 = 
Gc(^,t^'t^,t) — 1/2 then a trapezoidal approximation 
of the integrals results. If we call the matrix Gc, with 
1/2 subtracted from its diagonal elements, as Gc (note 
that we can treat one of the three independent momenta 
involved in Xc as a variable Q outside the matrix struc- 
ture), then we can write the two-particle Green function 
in a matrix form 



(Fq=o<^ctFq=o 



(60) 



(FqGc^Fq 



^^g-i(K,-Q).X,- 



where we have 



where (Fq)^ 

chosen i and j to index the cluster momentum-frequency 
space. This measurement may be performed efficiently if 
the product of three matrices in each set of parenthesis is 
tabulated as two sequential matrix-matrix products and 
stored before the direct product between the terms in 
parenthesis is calculated. When done this way, the cal- 
culation time required for this process scales like (NcNi)^ 
rather than {NcNi)'^ as would result from a straight- 
forward evaluation of the sums implicit in Eq. |6C| . 

For the reasons discussed above, Eq. |6^ becomes un- 
reliable at high frequencies \uJn\ > tt/At. The high fre- 
quency behavior of the two particle Green function can 
be recovered by using a method similar to that developed 



] . The first term of 
e diagram, is used 



for the one particle Green function |34 
its perturbation expansion, the bubb 
for the conditioning. It is calculated in two ways: First 
it is formed from the square of the properly conditioned 
cluster Green function. Second, it is calculated using the 
same approximation to the Fourier transform employed 
in Eq. The difference of the two may be used to 
condition the estimate 

Xc,,(Q) = 



Fq^o^^ciFq^o) (FqGcxFq 



FT _n ( G 



+Gc{K,)G,*iK^-Q)S,^ 



'Q 



(61) 



Moreover, this appends the right asymptotic behavior of 
the perturbation result to the two-particle Green function 
at high frequencies where QMC results are dominated by 
statistical errors. 
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E. Optimizing the Code 

In this subsection, we will discuss the optimization and 
parallelization of the QMC code. 

We generally find that the heat bath algorithm is more 
efficient, presumably because it has a lower acceptance 
rate and therefore deemphasizes the expensive step of 
updating the Green function. 

We may greatly reduce the statistical error in many of 
the measured Green functions by employing the trans- 
lational and point-group symmetries of the cluster. The 
QMC averaging over the HHS fields systematically re- 
stores the translational invariance of the system in time 
and space. So we may reduce the statistical error in the 
measured Green functions by averaging over all equiv- 
alent differences in spatial and temporal cluster coordi- 
nates. To reduce the statistical error further, we then 
average over all the lattice point group operations. For 
example, for Gc(K,a'„) 

Gc(K, iuJn) = y] Gc {TZ{K), tLOn) (62) 

where TZ is one of the symmetry operations in the point 
group of the lattice and N-ji is the total number of such 
symmetry operations. 

The two-particle Green functions typically have more 
statistical noise than their single-particle counterparts, 
and their matrices can be quite large. To reduce both the 
storage needed for these measurements and their statis- 
tical noise, the point group symmetry of the lattice may 
again be used. We first average the two-particle cluster 
Green functions over the different point-group operations 



Xo 



,(<7,K,K')--^^5]x.,.'(g,7^(K),7^(K')) . (63) 



We should also average over the symmetries of the di- 
agrams. I.e., for the particle-particle channel there are 
additional symmetries of the diagrams which include hor- 
izontal (K, iLj„; K', icjn') — > (K', ia;„'; K, ia;„) and ver- 
tical (K, zcj„; K', iwn') (— K, — ia;„; — K', — tj„') reflec- 
tions. After these symmetries have been imposed, we will 
lose no information and significantly reduce the storage 
requirements if we store x<t,(t' (q, K, K') for either K or 
K' within the irreducible wedge (we may not take both 
K and K' within the irreducible wedge though). 

The memory required for these calculations may be fur- 
ther minimized by limiting the use of the double precision 
arithmetic. The Green functions and all of the equa- 
tions associated with the calculation of the initial Green 
function and the Green function update, Eq. are com- 
puted with double precision (8 byte real) to minimize the 
problems with the accumulation of numerical error dis- 
cussed in Sec. [V B. However, to save memory, it is conve- 



nient to both calculate and store the two-particle cluster 
Green functions with single precision (8 byte complex). 
Since these measurements typically have a fraction of a 



percent statistical error, higher precision arithmetic and 
storage will not improve the accuracy of the two-particle 
measurements. 

The required CPU time may be reduced by optimiz- 
ing the inner loops. The two numerically most expen- 
sive parts of the QMC code are the Green function up- 
date, Eq. and the two-particle measurements, Eq. |6^. 
These can be written in terms of highly-optimized BLAS 
calls_[|D, DGER and CGEMM, respectively. To see that 
Eq. |49| can be calculated with an outer product, we define 



1) 



l + (l-G,,„,„0(e-"-(«"-«'")^l) 



Then Eq. k3 takes the form 



Gr n 



caij 



~l~ i^caim ^i,m)o,YnG car. 



(64) 



(65) 



This is a vector outer product and matrix update, with 
vectors a„i{Gcaim — ^i,m) and Gcamj for fixed to. 

Additional speedup of the calculation is possible by 
writing parallel codes. The DCA-QMC codes are ex- 
tremely well suited for massively parallel supercomputers 
because of their efficient use of the floating-point capa- 
bilities of such machines and the highly parallel nature of 
the codes and the underlying algorithm. With the cur- 
rent relative decline in the availability of vector super- 
computers and the concomitant increase in the number 
of massively parallel supercomputers, this is an impor- 
tant feature of the algorithm. In the remainder of this 
subsection, we discuss flrst the general parallel nature of 
the algorithm. 

There is a high degree of parallelism in the DCA-QMC 
algorithm, which one may exploit. This parallelism exists 
on two levels. First, QMC is itself inherently parallel 
because it consists of a number of stochastic random- 
walks. One may think of QMC as one long Markov- 
chain walk. Measurements are made periodically along 
this walk. At the end of the walk, these measurements are 
averaged and the flnal result, with error bars, is obtained. 

However, there is no reason why this Markov-chain 
walk has to be continuous. It has been known for 
years that one can perform several independent, shorter 
Markov-chain walks and average the results of each walk 
to obtain the flnal result of the calculation. The result 
can be an almost perfect parallel speedup as an increasing 
number of processors is applied to a problem. This arises 
because only an extremely small amount of communica- 
tion between processors is required - flrst to initialize 
the Markov-chain walks and then to collect the data for 
averaging at the end of the Markov process (even this 
averaging can be done in parallel using MPI calls). We 
call this the "Perfectly Parallel" algorithm. 

The second degree of parallelism exists in the linear al- 
gebra problem itself. That is, one can distribute the vec- 
tors and matrices which comprise the linear algebra prob- 
lem across several processors. (The matrix in our case is 
the Green function discussed above.) Such a break-up 
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of the data becomes of paramount importance when the 
size of a matrix is so large that it cannot possibly fit 
within all of the memory available on a single processor 
of a computer. 

The issue of interprocessor communication now be- 
comes paramount as one performs linear algebra. How- 
ever, two things work in our favor here. First, the main 
linear algebra operation of the QMC is a vector outer 
product - which is in itself inherently parallel. Second, 
this is a well-studied problem and again an efficient li- 
brary package, the parallel PBLAS Q, exists to solve 
it. When we divide the Green function over all of the 
processors that we use in a run on a parallel machine, we 
call this the "Truly Parallel" method. 

The Truly Parallel method can be used to efficiently 
fill all available processors of a parallel machine with one 
DCA-QMC problem. Often, however, the problem of 
interest is not so big as to require the entire machine 
for one Green function, but is too big to fit within the 
RAM available on a single processor and hence too big 
for the Perfectly Parallel code. To efficiently use available 
hardware for these problems, one can employ a "Hybrid" 
code, which is both truly parallel in part and perfectly 
parallel in part. 

The hybrid code may be thought of as using blocks of 
processors to distribute Green functions and in turn per- 
forming a perfectly parallel QMC with many such blocks. 
For example, say that the Green function for the prob- 
lem at hand will not fit in the memory of a single pro- 
cessor, but will fit within the memory of 4 processors. 
Assume also that there are 100 processors available for 
a run. The hybrid code then allocates all 100 proces- 
sors, divides these 100 processors into 25 blocks of 4, 
distributes copies of the initial Green function onto each 
of the 25 blocks, and then does a perfectly parallel QMC 
using these 25 blocks. This makes the most use of the 
resources of a machine and is especially well-suited for a 
Symmetric Multi-Processor machine, where many nodes 
exist and each node comprises several processors with a 
shared, relatively large, pool of RAM. 



V. THE DCA ALGORITHM 

In this section, we will discuss how the QMC and DCA 
formalism are combined into a DCA algorithm for simu- 
lating lattice problems. The complete DCA program is 
made of three completely separate parts as illustrated in 
Fig. ^ The first part is the self-consistent loop which is 
the main part of the algorithm. It includes the DCA self- 
consistency loop composed of the QMC block and the 
coarse-graining of the lattice. This program is usually 
run on a parallel supercomputer, and the cluster self en- 
ergy and the various two-particle cluster Green functions 
are written into files. In the second part various one- 
particle and two-particle lattice Green functions are cal- 
culated from the cluster Green functions obtained from 



the self-consistent loop. This part of the code is generally 
run on a workstation and it requires in the data gener- 
ated by the first part of the code. The third is devoted to 
the analytical continuation of the imaginary-time Green 
functions to real frequencies using the Maximum Entropy 
Method (MEM). 



A. Part 1: The self-consistent loop 

1. The DCA iteration procedure is started by setting 
the initial self-energy Sc(K,iw„) — 0, or to a per- 
turbation theory result. 

2. E is then used to compute the coarse-grained Green 
function G(K,iw„), 

G{K,zLo^)^Y.-- ^ V . w (66) 



3. The next step of the iteration is to use G'(K,ia;„) 
to compute the host Green function Q{K, iuJn)^^ — 
G(K.,i(jJn)~^ + 5]c(K, zcj„) which must be intro- 
duced to avoid over-counting diagrams. C?(K,iaj„) 
serves as the input to the QMC simulation to yield 
a new estimate for the cluster self-energy. 

4. t/(K, iLUn) rnust be Fourier transformed from 
the momentum-frequency variables to space- 
imaginary-time variables before being introduced 
in the QMC part of the program as the initial Green 



function Gc 



^3 J 



corresponding 



to all Si — 0. Eq. is used to generate the cluster 
Green function corresponding to = 1, or to the 
{si} from a previous run. 

5. The QMC step is next and is the most time con- 
suming part of the algorithm. Each QMC step is 
warmed up before one starts to perform measure- 
ments. While making measurement, we average 
over the differences in space and time and the point 
group operations, as described above, to reduce 
the statistical error. This together with the QMC 
averaging restores the translational invariance of 
the system in time and space, so {Gcij}^^.-^ = 

Gc(Xi — X.j,Ti — Tj). 

6. G'c(Xi — XjjTi — Tj) is then Fourier-transformed 
to GcCKjiujn)- We calculate a new estimate 
for the self-energy Ec(K, icjn) = Q(K,iuJn)~^ — 
GciK,iuj^)-\ 

7. Starting with step 2, the procedure is repeated until 
I]c(K, icjn) converges. This typically happens in 
less than ten iterations. The number of iterations 
decreases when Nc increases since the coupling to 
the host is smaller {0{1/Nc)) § for larger clusters. 
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The convergence test is made on the ratio p, 

I EK(^c„«„(K,iwo) - Sco;d(K,«Wo))| 



IEK^co/d(K,ZCJo)| 



(67) 



where ujq — nT 



8. Once convergence is reached to the desired accu- 
racy, the remaining one and two-particle measure- 
ments are made in a final QMC iteration. As in 
a usual QMC simulation, bins of measurements 
are accumulated and error estimates are made 
from the fluctuations of the binned measurements. 
These error estimates are accurate provided that 
the bins contain enough measurements so that the 
bin-averages are uncorrelated. The statistical er- 
ror may be reduced by averaging over the different 
symmetries as discussed above. 

Once the cluster Green functions are obtained, the de- 
termination of the lattice quantities requires additional 
steps which are done in separate programs. 



QMC 

G.(x), x^(T) 



g'^= G" + E 



G(K) = I G(K + k ) 
k 



MEM 
N(co), x(co) 



Analysis Code 
Z(T), n(k) 



FIG. 7. Sftctcli of tfio DCA algorithm: the self-consistent 
loop, the analysis part and the MEM part. 



positive-definite function remains everywhere positive). 
We also use a multidimensional spline interpolant, like 
some Akima splines, which does not overshoot. However, 
it is important to note that this interpolated self-energy 
should not be used in the self-consistent loop as this can 
lead to violation of causality 

The interpolated self energy, S(k, iujn) is then used to 
calculate the Fermi surface. For this we use the discrete 
form of Vn(k) 



A«k _ T E„ G(k + Ak, tLJn) - G(k, tLJn) 



Ak 



Ak 



(68) 



which, in a Fermi liquid (or a marginal Fermi liquid) is 
maximum at the Fermi surface. The quasiparticle weight 
may be approximated with 



1 



(69) 



which becomes exact as T 0. 

The calculation of the lattice susceptibilities in the 
particle-hole channel and in the particle-particle chan- 
nel is also made in this code. The stored QMC cluster 
susc eptibi lities are used for this purpose as prescribed in 
Sec. [HE . Here, we first form the corresponding coarse- 



grained and bare cluster susceptibilities, and then we use 
Eq. ^ to calculate the corresponding coarse-grained lat- 
tice susceptibility. To calculate the pair field susceptibil- 
ities, we first calculate the corresponding coarse-grained 
two-particle reducible vertex, and then use Eq. ^ to cal- 
culate the coarse-grained lattice pair-field susceptibility 
matrix. 

In the two-particle calculations, it is tempting to in- 
terpolate the cluster vertex functions to the lattice mo- 
menta. However, this would increase the size of the ma- 
trices which must be inverted in Eq. ^ dramatically, 
making the calculation of the lattice susceptibilities nu- 
merically much more expensive. 



B. Part 2: Numerical calculation of lattice quantities 

The self-consistent loop yields cluster Green func- 
tions Gc(K, zcij„), I]c(K, icjn), and susceptibilities 
Xco-,o-'(K,iWri;K',za;m), Hg(K, iw„; K', zw^) which may 
be used to construct the equivalent lattice quantities. 
This is done in a separate computer program in which the 
irreducible quantities of the cluster which are in the DCA 
approximation identical to those of the lattice are used to 
compute the corresponding reducible lattice quantities. 

To calculate the single-particle quantities, an interpo- 
lated self-energy i;(k, zcj„) may be used. This is espe- 
cially important for the calculation of |Vn(k)|, and other 
quantities such as band structure, where continuity of the 
self energy is important. We often use bilinear interpo- 
lation for this purpose, since it is guaranteed to preserve 
the sign of function (i.e. the bilinear interpolation of a 



C. Part 3: Analytic continuation 

Unfortunately, there is no reliable way to perform the 
direct analytic continuation of Sc(K, iuJn)- Fade approx- 
imants lead to very unstable spectra because of the QMC 
statistical noise contained in Ec(K, iwn). The binned 
imaginary-time Green function data accumulated from 
the cluster calculation must be used to obtain lattice 
spectra from which Ec(K,w) may be deduced. To ob- 
tain the self-energy and spectral-weight function A(k, lu) 
of the lattice in real frequencies, we first compute the 
cluster spectral- weight A(K,t£). This is done using the 
Maximum entropy method p5[ to invert the following 
integral equation 

GiK,T)^ J dio^^^A(K,Lu) , (70) 
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where G(K, r) is the imaginary-time Green function ob- 
tained from the QMC simulation of the cluster. 

Since j4(K,a;) ~ — l/7r/mG'(K, w), the full frequency- 
dependent coarse-grained Green function G(K,uj) is ob- 
tained using Kramers-Kronig relations. Then, the equa- 
tion 



G(K,c.) 



N 



E 



1 



uj — e — e 



K+k 



(71) 



is solved for the real-frequency self-energy Sc(K, lu) using 
a complex root finder [3^ . This self energy may then be 
inerpolated onto the lattice k points using a high-level 
interpolant which also preserves the sign of the imaginary 
part. 

The above steps are unnecessary if the local quantities 
are to be computed since the local lattice and cluster 
Green functions correspond one-to-one. For example, we 
may directly analytically continue the local cluster Green 
function to obtain the lattice density of states. 



VI. APPLICATION TO THE 2D HUBBARD 
MODEL 

We will now show the results of the application of the 
DCA to the two-dimensional Hubbard model. The Hub- 
bard model has a long history and is believed to con- 
tain the mechanism of various physical phenomena such 
as magnetism, metal-insulator transitions and more re- 
cently superconductivity and non-Fermi liquid behavior. 
Our intent in this section is not to exhaustively study this 
model's properties, but rather to use it to illustrate the 
power and limitations of the DCA and to survey what 
can be done. 

Since the two-dimensional model is not expected to 
have a finite-temperature magnetic or perhaps even su- 
perconducting transition, we will add a hopping t± |3^] 
into the third dimension between an infinite set of weakly 
coupled Hubbard planes 

t±{kx, ky,kz) = — 2tj_(cos kx — cos kyf' cos k^ ■ (72) 

We take t± <C t, and treat the additional coupling at the 
DMFA level, so the self-energy is independent of kz- This 
is accomplished by modifying the coarse-graining cells 
into rectangular solids of dimensions Afc, Afc and 2t: in 
the kx, ky and k^ directions, respectively. After coarse- 
graining, the problem is reduced to a two-dimensional 
cluster. Information relevant to the mean-field coupling 
between the planes is contained within Q. 



A. Results at Half-filling 

The physics of the half filled model is a severe test of 
the DCA as well as finite-sized simulations (FSS) due to 
the quantum critical point at zero doping. As this point 



is approached, both the dynamical and spatial correlation 
lengths diverge, and both the DCA and FSS are expected 
to fail. 



1. Antiferromagnetism 

Earlier finite size simulations ||3^ , |39| employing the 
QMC method have led to the conclusion that the ground 
state is an antiferromagnetic insulator at half-filling. 
Since the model is two dimensional, we know from the 
Mermin- Wagner theorem that the transition temperature 
is necessarily zero. But as found in infinite dimensions 
p2t , the DMFA predicts a finite temperature transition 
even in two dimensions. This spurious behavior may be 
attributed to the lack of non-local correlations in the 
DMFA. These correlations are known to induce strong 
fiuctuations particularly in reduced dimensions and are 
responsible for the suppression of the finite temperature 
transition. The DCA which includes these non-local cor- 
relations is thus expected to progressively drive the spu- 
rious finite temperature transition found in the DMFA 
towards zero temperature as the cluster size increases. 



0.5 



0.4 



, 0.3 



0.2 



0.1 



0.1 
0.09 
0.08 
0.07 
-0.06 
0.05, 











"o 

, 1 , 





10 



0.05 




1.5 N=1.00 






N 


=1 T^=0.094-)fcl.04 


□ 


N 


=8T^=0.087-)fcl.06 


o 


N 


= 16 T =0.078 7=1.08 

N • 


A 


N 


= 18 T =0.074 7=1.12 

N • - 


* 


N 


=32X^=0.073 7=1.11 


X 


N 


=36X^=0.068 7=1.14 " 



FIG. 



0.15 

T 



0.25 



8. The inverse antiferromagnetic susceptibility versus 
temperature for various cluster sizes. The lines are fits to the 
function (T — Tn)'*'. In the inset, the corresponding Neel tem- 
peratures, determined by the divergence of the susceptibility, 
are plotted. The line is a polynomial fit to the data, excluding 
Afc = 4. 

This behavior is illustrated in Fig. where the in- 
verse antiferromagnetic susceptibility is plotted versus 
temperature for (5 = and various values of Nc which 



preserve the lattice symmetries as discussed in Sec. Ill C . 
At high temperatures, the susceptibility is independent 
of Nc, due to the lack of non-local correlations. In con- 
trast to FSS calculations, the low temperature suscepti- 
bility diverges at T = Tn, indicating an instability to an 
antiferromagnetic phase. As Nc increases Nc > I, the 
non-local dynamical fluctuations included in the DCA 
suppress the antiferromagnetism. For example, when 
Nc = 1, the susceptibility diverges with an exponent 
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7 ~ 1, as expected for a mean- field theory; whereas the 
susceptibilities for larger values diverge at lower tem- 
peratures with larger exponents indicative of fluctuation 
effects ^l|. At first, these effects are pronounced; how- 
ever, as Nc increases, Tn falls and 7 rises more slowly 
with increasing N^. This can be understood from the 
singular nature of the spin correlation length, which at 
least in the large U limit is expected to vary as ^ oc C^/"^, 
where A is a constant of the order of the exchange cou- 
pling. For this quantum critical transition, we expect the 
DCA to indicate a finite temperature transition once ^ 
exceeds the linear cluster size. Since correlations build 
exponentially, large increases in the cluster size will only 
reduce the DCA transition temperature logarithmically. 

Note that the data for Tn(A^c) falls on a smooth curve, 
except for T^{Nc = 4). This behavior was seen pre- 
viously in the transition temperature of the Falicov- 
KimbaU model, calculated with DCA. The ^ 4 
data falls well off the curve produced by the other data, 
and has a much larger exponent indicating that fluctua- 
tion effects are more pronounced. Presently this behavior 
is not completely understood but may be related to the 
fact that the maximum coordination number for Nc = 4 
is two, whereas it is greater than two for cluster sizes 
larger than iVc = 4. 
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FIG. 9. Tn versus Nc for different values of t±/t 

An interplanar coupling can significantly alter the 
phase diagram. However, since the superexchange cou- 
pling varies roughly like the square of the hopping, it is 
necessary to make tj_ a significant fraction of the intra- 
planar coupling t in order to see an effect. For example, 
if t±/t = 0.4, the ratio of the interplanar to intraplanar 
exchanges is roug hly J^/J w 0.16. In Fig. | the antifer- 
romagnetic transition temperature is plotted versus Nc 
when tA^/t = 0,0.4,1.0 when U = W = 2 and 6 = 0. 
For both t±/t — 0.4 and t±/t — 1.0, the transition tem- 
peratures for A'c = 16 and 32 are the same to within the 
numerical error. Thus, the finite-temperature transitions 
found for small clusters, can be preserved as Nc — > 00 by 
introducing the interplanar coupling. 



2. Mott transition at half-filling 

In the strong coupling limit, a Mott Hubbard gap is 
expected to open in the charge excitation spectra. In the 
weak coupling limit, the situation is less clear. Since the 
ground state of the half filled model is always an antifer- 
romagnet, the system remains insulating, but the nature 
of the insulating state in weak coupling is less clear, and 
depends upon the dimensionality. In one dimension, Licb 
and Wu |^ showed long ago that a charge gap opens as 
soon 'AS U > 0. There is a spin-charge separation and 
there is no long range order, even at T = 0. Hence, the 
Slater scenario is not responsible for the metal-insulator 
transition and the low energy spin excitations are de- 
scribed by the Heisenberg model. In infinite dimension, 
the model can be mapped to a self-consistent Anderson 
impurity problem. The solution of the self-consistent 
equations have been obtained numerically by various au- 
thors. For small U <C W, the antiferromagnetic transi- 
tion temperature Tn is higher than any temperature at 
which there is a metal-insulator transition in the param- 
agnetic phase. Hence, the metal-insulator transition in 
infinite dimension is due to the Slater mechanism. In 
two dimensions, the Mermin- Wagner theorem prohibits 
long range order for any T > 0, and we find that the 
weak coupling transition is similar to what is found in 
one dimension. 
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FIG. 10. The single-particle density of states p{u) at 
/3 = 32, [/ = 1, and t± = 0. 



The density of states p{uj) (shown in Fig. 10) confirms 
the destruction of the Fermi liquid quasi-particle peak 
by short-range antiferromagnetic correlations. With in- 
creasing Nc, the gap opens fully, and the Hubbard side 
bands become more pronounced. 

In Fig. ^ the behavior of Tn is compared to the tem- 
perature Tg where the gap opens. In contrast to Tn, 
Tg increases with the size of the cluster. This confirms 
the conclusion that a gap, which is not due to antiferro- 
magnetism alone, opens at finite temperatures in the 2D 
Hubbard model. 
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FIG. 11. The Neel, Tn, and gap, Tg, temperatures versus 
Nc when U = W = 2, S = and tj_ = 0. 

3. Comparisons with finite system simulations 

A great number of simulations have been performed in 
the half-fiUing regime of the two-dimensional Hubbard 
model over the last decade. Most of them are based on 
QMC with imaginary-time data analytically continued 
by the maximum entropy method. While these studies 
all agree for U > W, for U <W they have led to conflict- 
ing results The reason is that the metal-insulator 
transition is related to the antiferromagnetic transition so 
that it is difficult to distinguish between the two physical 
processes. As a consequence various conflicting scenar- 
ios for the disappearance of the quasiparticle peak at low 
temperatures have been proposed. These controversies 
are inherent to the limitations of the finite system simu- 
lations. There are artificially introduced finite size gaps 
when the correlations become comparable to the system 
size. This does not occur when U > W because the Mott 
gap opens well before the magnetic correlations set in. It 
is thus fair to ask to what extent the conclusion reached 
with the DCA above may be more reliable. For this it is 
necessary to compare the DCA to FSS. 

In Fig. |l^ and Fig. |l^ we show the imaginary-time 
Green function G{t) at the Fermi point X = (tt, 0). This 
quantity has a more rapid decay from its maximum at 
G{P/2) when the effects of the correlations are stronger. 
In finite systems, the decay is sharper for smaller lat- 
tices while in the DCA it is the opposite. This behavior 
marks the fundamental difference between the FSS and 
the DCA. At low temperatures, in FSS, the correlation 
length is greater than the lattice size. Thus, the effects 
of the correlations are overestimated for smaller clusters 
because these systems are artificially closer to criticality 
than a system in the thermodynamic limit. This ten- 
dency is reduced by increasing the cluster size, which 
moves the system in the direction of the thermodynamic 
limit. The situation is radically different in the DCA 
where the system is already in the thermodynamic limit. 
The DCA approximation consist in restricting correla- 



tions to within the cluster length in the infinite system. 
As the cluster size increases, possible longer-range cor- 
relations are progressively included. Thus, the effects of 
the correlations increase with the cluster size. 
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FIG. 12. The imaginary-time Green function at the point 
X = (7r,0) on the Fermi surface from finite size QMC (filled 
symbols) and from DCA (open symbols) with U = 1.1, 
/3 = 16, tx = and Nc — 16, 36, 64. The size increases 
from top to bottom for FSS. The size increases from bottom 
to top for DCA. 

The density of states shown in Fig. |l^ supports these 
conclusions. The finite size gap in the FSS decreases 
when the cluster size goes from Nc = 16 to 64. While, 
in the DCA for Nc — 16 there is a pseudogap that turns 
into a true gap when the cluster size is increased to 64. 
Since by construction, the DCA underestimates the gap, 
we can affirm that at this temperature, the gap exists in 
the thermodynamic limit. Its actual value is bracketed 
by the FSS and the DCA. This behavior is characteristic 
of the DCA. It has been extensively verified on the one- 
dimensional Hubbard model Q . 



FSS Ai 

-f!^ \ ^ 


ll 
■J 

V .1 J 


N=16 

N =64 


DCA /, 





-2-10 1 2 



CO 

FIG. 13. The density of states p{uj) from finite size QMC 
(top) and from DCA (bottom) at [/ = 1, /3 = 32, = and 
Nc = 16,64. 
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B. Results away from half- filling 

1. Sign Problem 
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FIG. 14. The average sign as function of the inverse tem- 
perature /3 for iVc = 8 at 5 = 0.1 for U = 1.0, 1.5,2.0. In 
the inset, the average sign is plotted versus doping 6 when 
(7 = = 2, = and /3 = 54. 

The most serious limitation of QMC calculations at 
low temperatures is the sign problem. Off half-filling, 
the sign of P({s,}) oc det[Gl{{s,})] x det[G^^{{s^})] can 
be negative so that it can no longer be interpreted as a 
probability distribution. The solution is to reinterpret 
\P{{si})\ as the probability of the configuration {si} and 
associate its sign with the measurement. For any 
operator O, this becomes 



(O) 



sign({si}) 



(73) 



where sign({si}) is the sign of P{{si}), 0{{si}) is the 
value of the operator for the field configuration {s^}, and 
the primed sums are over configurations generated by im- 
portance sampling. In finite system simulations, as the 
temperature is lowered, the average sign becomes expo- 
nentially small 1^,0 so that it is no longer possible to 
obtain good statistics. This sign problem has posed a 
formidable challenge in the field of numerical simulations 
for nearly two decades. 

Some recent works have brought some hope. Guber- 
natis and Zhang |43| and Zhang |44| have shown that by 
putting a constraint on the fermion determinant, one can 
construct an approximate algorithm which shows some 
improvement on this problem. While the resulting algo- 
rithm seems to be free from the sign problem, it is possi- 
ble that the constraint introduced may affect the ergodic- 
ity of the algorithm. The ergodicity question is suggested 
by the work of Sorella who employ a similar idea as 
the former authors but who arrived at different results. 
The most promising new direction seems to be that of 
Chandrasekharan and Weise ]46| . They proposed a new 
algorithm which is rigorously free from a sign problem for 
certain classes of models. The basic idea is to decompose 



a configuration of fermion world-lines into clusters that 
contribute independently to the sign. There are two type 
of clusters: clusters whose flip changes the sign called 
meron and others that do not modify the sign after a flip. 
Configurations containing mcron-clusters contribute to 
the partition function, while all other configurations con- 
tribute 1. Hence, this cluster representation describes the 
partition function as a gas of clusters in the zero-meron 
sector. 

The sign problem remains in DCA simulations, as il- 
lustrated in Fig. |l^ where the average sign is plotted 
versus inverse temperature for various values of U when 
S = 0.1 and Nc = 8. In the inset, the average sign is 
plotted versus doping when U = W = 2, jS = 54 and 
Nc — 8. As in FSS, the sign is worst just off half fill- 
ing. However, the DCA sign problem is significantly less 
severe than that encountered in FSS. This is illustrated 
in Fig. |l^ where the average sign for the DCA and the 
BSS simulations of White et al. are compared when 
U/W = 1/2, (5 = 0.2, and tx_ = 0. When tx_ is finite, the 
average sign increases further (not illustrated). We at- 
tribute this strange behavior to the action of the host on 
the cluster, but its actual mathematical justification is 
still mysterious. Nevertheless, due to the large reduction 
in the severity of the sign problem, we are able to study 
the physics at significantly lower temperatures than is 
possible with FSS! 
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FIG. 15. A comparison of the average sign for the DCA 
and FS simulations H when U/W = 1/2,5 = 0.2, ti_ = 0. 
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2. Single-Particle Properties 

Much can be learned about the single-particle prop- 
erties of the system, especially Fermi liquid formation, 
from studying the momentum distribution function n(k), 
the single-particle spectra A(k, w) and the single-particle 
self energy E(k, w). For a Fermi liquid, the self energy 
S(kF,w) ~ (1 - l/Z)Lv-ibuj'^ where 6 > 0, 1/Z > 1 and 
ki? is a point on the Fermi surface. The corresponding 
^(ki?,^) is expected to display a sharp Lorentzian-like 
peak, and |Vri(k)| is also expected to become sharply 
peaked at the Fermi surface. In each case, these quan- 
tities are calculated by first interpolating the cluster self 
energy onto the lattice k points. 

For example, the gradient of the momentum distribu- 
tion function is plotted in Fig. |l6| when U — 1, (3 = 44, 
S = 0.05 for different values of Nc (this temperature 
would correspond to roughly room temperature for the 
cuprates in units where the bare bandwidth W — 2eV). 
Apparently, at this temperature, there are two Fermi sur- 
face features, one centered at F = (0, 0) and one centered 
at M = (7r,7r). The Fermi surface centered at F = (0,0) 
has roughly the volume expected of non-interacting elec- 
trons, so we will call it the electron-like surface and the 
other hole-like. Note that the hole-like Fermi surface be- 
comes more prevalent, and the peak near (7r/2,7r/2) di- 
minishes, as Nc increases. We therefore attribute this 
behavior to short-ranged correlations. 
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X 

FIG. 16. |Vn(k)| versus k when [/ = 1 /3 = 44, fx = and 
S = 0.05 for = 1, 8 and 16. 

We can further resolve the different surface features, 
by investigating the single-particle spectrum A(k, uj) as 
shown in Fig. |l^ for f7 = 1, /5 = 44, (5 = 0.05 and 
Nc — 16. The graph (d) on the upper left of Fig. |l^ plots 
the corresponding location of the maxima of |Vn(k)|. 
Along the direction from F to ill, A(k, oj) shows a rela- 
tively well defined and symmetric peak at lu — ai the 
location as indicated by the maximum of |Vn(k)|. The 
only notable feature is that the peak is a bit better de- 
fined for k closer to the zone center F. Along the direction 
from M to X, the part of the hole-like Fermi surface clos- 
est to the X point is resolved. Here the peak in A{'k,uj) 
crosses the Fermi surface at roughly the same k where the 
peak in |Vn(k)| is seen; however, the peak in A(k, w) is 
broader and is heavily skewed to higher frequencies. Fi- 
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nally, along the direction from X to F, we find very sharp 
peaks in ^(k, a;); however, none occur at a; = indicat- 
ing that there are well defined quasiparticle excitations 
along this direction, with a small pseudogap, presumably 
due to the short-ranged order. This pseudogap behavior 
becomes more pronounced as the temperature is lowered. 



=0.05 U: 




Q) CO 

FIG. 17. (a)-(c) The single-particle spectrum j4(k,tj) for 
f/ = 1, /3 = 44, (5 = 0.05, tj_ = and iV^ = 16 along cer- 
tain high symmetry directions. The arrows and bold lines in 
Figs, (a) and (c) indicate the spectra which cross the Fermi 
energy with a peak closest to w = 0. In (b), no such peak is 
found which crosses the Fermi energy, (d) the maxima of the 
|Vn(k)| data illustrated in Fig. ^ plotted versus k. 

This can be seen in the density of states, shown in 
Fig. ^ where the gap is more pronounced. At high tem- 
peratures, (3 — 4: the Hubbard side bands are apparent 
at w « ±1/2. As the temperature is lowered, a central 
peak begins to develop. At low temperatures, /3 > 24 a 
pseudogap begins to develop. 




FIG. 18. The single-particle density of states p{uj) when 
U = 1, S = 0.05, t± — and Nc = 16 for several different 
values of the inverse temperature 13. As the temperature is 
lowered, p{u!) develops a pseudogap due to the short-ranged 
antiferromagnetic order. 

More can be learned by investigating the self energy 
directly. In FigT9, both the real and imaginary parts 
of the self energy are plotted for the three values of k 
indicated by filled circles in Fig. [l7|(d). The self energy 
on the part of the Fermi surface along the direction from 
F to M, looks roughly Fermi-liquid like. However, the 
self energy on the parts of the Fermi surface closest to X 
have pronounced non-Fermi liquid character, especially 
at k = (tt, 0.48) where the real part displays a minimum 
and the imaginary part crosses the Fermi energy almost 
linearly. At k = (2.571, 0), the real part again displays a 
minimum, but the imaginary part has an almost Fcrmi- 
liquid-like maximum at the Fermi energy, and then once 
again the scattering rate increases dramatically at higher 
energies. All of the points close to X share this dramatic 
asymmetry; that excitations below the Fermi energy are 
much longer lived than those above. Thus, we expect 
that the transport from these parts of the Fermi surface 
would be predominantly hole-like. 
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FIG. 19. The low-frequency self energy, plotted versus uj 
for the three k-points denoted by filled circles in Fig. p^d) 
where the Fermi surface defined by the maxima of |Vn(k)| 
crosses the high symmetry directions. 

The non-Fermi liquid features, including the hole-like 
distortion of the Fermi surface, the anisotropy and non- 
Fermi liquid features of the self energy, and the pseudo- 
gap in the density of states, become more pronounced as 
Nc increases. Thus, it is reasonable to assume that these 
features will persist as A^c ~* oo. 



3. Superconductivity 

We searched for many different types of superconduc- 
tivity, including s, extended-s, p and d-wave, of both odd 
and even frequency and we looked for pairing at both 
the zone center and corner. Only the pairing channels 
with zero center of mass momentum (zone center) are en- 
hanced as the temperature falls. Of these, only the even- 
frequency d-wave pair field susceptibility diverges. This 
is illustrated in Fig. where all of the zone center sus- 
ceptibilities are plotted versus temperature for U = 1.5, 




FIG. 20. Various pair field susceptibilities calculated at 
the zone center and plotted versus temperature for U =1.5, 
8 = 0.05, t± = 0, and Nc = 8. Pairing is found only in the 
even-frequency q = d-wave channel. 

In contrast to the antiferromagnetic susceptibility 
which falls as Nc increases, the d-wave pair field sus- 
ceptibility generally rises with Nc, except at very low 
temperatures. This is illustrated in Fig. |^. However, 
for Nc = 16 at low T, it falls abruptly when T < 0.03. 
This behavior is consistent with the lack of superconduc- 
tivity in the purely two-dimensional model. However, 
in the inset, we see that a very small interplanar cou- 
pling t±/t — 0.2 causes the susceptibility to continue to 
rise with decreasing temperature. Thus, perhaps a very 
small interplanar coupling is able to stabilize the mean- 
field superconductivity seen in smaller clusters 1^,^ . 




FIG. 21. The d-wave pair field susceptibility versus T for 
several values of Nc when U = 1.5, S — 0.05 and t±/t — 0.0. 
In the inset the d-wave pair field susceptibility is plotted ver- 
sus T when Nc = 16, (7 = 1.5 and S = 0.05 for tj_/t = 0.2. 
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VII. SUMMARY 

We have presented the algorithmic details of the dy- 
namical cluster approximation. The technique consists 
of mapping the original lattice problem to a problem in- 
volving a finite cluster dynamically coupled to an infinite 
host. The cluster problem may be solved by a variety 
of techniques that include the QMC method, the FLEX 
approximation or the NCA. 

An extensive account of a QMC method used to solve 
the cluster problem was given. Though this algorithm 
requires significantly greater computer power than the 
Blanckenbecler-Sugar-Scalapino algorithm which is often 
used for finite systems, it has some advantages. First, it 
does not show any numerical instability at low tempera- 
ture; thus, it avoids the time costly matrix factorization 
step that slows down the BSS algorithm. Second, the al- 
gorithm is quite general and can be applied to problems 
that do not have any explicit Hamiltonian formulation 
with known parameters. Third, a mean-field coupling 
to other degrees of freedom may be easily incorporated. 
Fourth, the minus sign problem is far less severe in DCA 
simulations. This allows us to study systems at signifi- 
cantly lower temperatures, with stronger interactions, or 
with larger clusters than can be studied with the BSS 
algorithm when the sign problem is apparent. 

The full DCA algorithm is made of three separate units. 
In the first unit the coarse-graining of the lattice is per- 
formed and the resulting self-consistent cluster problem 
is solved by the QMC technique. This unit requires the 
formidable computer power available on massively par- 
allel computers. The second part deals with the calcula- 
tion of the lattice one and two-particle Green functions 
from those of the embedded cluster. In the last part 
the analytical continuation of the imaginary-time Green 
functions to real frequency Green functions is performed. 

In order to illustrate the originality of the DCA tech- 
nique, we have applied it to the two-dimensional Hubbard 
model with a small interplanar mean-field coupling. 

The DCA method was developed to address some of the 
shortcomings encountered in the dynamical mean-field 
theory. The lack of non-local fluctuations in the DMFA 
leads to incorrect predictions when this method is applied 
to systems in finite dimension. In particular, we have 
seen that in violation of the Mermin- Wagner theorem, 
the DMFA predicts a finite-temperature transition in the 
two-dimensional Hubbard model. We have shown that 
in the DCA, this transition is progressively suppressed 
as the range of the fluctuations (i.e the cluster size) is 
increased. 

We also find that a finite-temperature gap persists well 
into the weak coupling regime of the half filled model. 
As the DCA systematically underestimates the gap for- 
mation, these conclusions are valid in the limit Nc — *■ oo. 
Since the temperature where the gap opens increases with 
Nc, while Tn decreases, the Slater mechanism is likely not 
responsible for the metal-insulator transition in the two 



dimensional Hubbard model. The resulting phase dia- 
gram is consistent with Anderson's view that the effective 
Hamiltonian for the 2D Hubbard model at half-filling for 
all [/ > and A ^ T (where A is the gap energy) is the 
2D Heisenberg Hamiltonian p9| . 

We find no evidence for a Kondo peak, or the associated 
Fermi liquid behavior for the unfrustrated model near 
half filling. Since this is an essential feature of the DMFA 
solution of the doped model or the half filled model with 
U ^ W, we conclude that the DMFA is a very poor 
approximation for the two-dimensional model, especially 
for behavior such as the Mott transition, observed near 
or at half filling. 

When the model is doped, the sign problem becomes 
significant and will certainly affect the quality of the re- 
sults. However, the sign problem is significantly less se- 
vere than that found in finite size systems, allowing us to 
explore these model systems at significantly lower tem- 
peratures, larger coupling or larger clusters than hereto- 
fore possible. In the doped model, we find evidence of 
non-Fermi liquid behavior even for relatively small val- 
ues of U/W. This has been observed in the self energy, 
single particle spectra, and density of states. We also 
find that the d-wave pair-field susceptibility is divergent 
for small clusters. A trend that is not present in the 
DMFA because the method cannot treat non-local order 
parameters. 

Finally, the DCA is a very versatile technique that may 
be applied to a variety of problems. A straightforward 
generalization of this algorithm to the periodic Ander- 
son model in two dimensions will allow us to study the 
physics of the recently discovered two-dimensional heavy 
fermion systems. In its present form, we have incorpo- 
rated diagonal disorder in the 2D Hubbard model which 
will allow us to address the interesting problem of disor- 
der and interaction in two dimensions. A future improve- 
ment of the DCA algorithm itself is to insert long range 
fluctuations in the algorithm which would be treated per- 
turbatively. 
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